Midway Experiment

Prepare the data

Load the packages

library(lme4)
library(lmerTest)
library(afex)
library(effects)
library(car)
library(MuMIn)
library (dplyr)
library(emmeans)
library(DHARMa)
library(performance)
library(patchwork)
library(rstatix)
#for plots and tables
library(ggplot2)
library(ggpubr)
library(forcats)
library(RColorBrewer)
library(tidyverse)
library(sjPlot)
library(sjmisc)
library(gt)
library(purrr)
library(stringr)
library(tidyr)
library(hms)
library(lubridate)
library(ggeffects)

Load this file for normalized to quantum efficiency of photosynthesis per Silsbe and Kromkamp (2012)

mid_ps <- read.csv("../../../data/midway_2024/transformed/midway_ek_alpha_normalized_2024.csv")

Rename specimenID to simply ID

mid_ps <- mid_ps %>% 
  rename(ID = specimenID)

Standardize date to ISO

mid_ps$date <- ymd(mid_ps$Date)

Make sure time is hms

mid_ps$rlc_end_time <- as_hms(mid_ps$rlc_end_time)

Assign run as a factor

mid_ps$run <- as.factor(mid_ps$run)

Assign temperature as a factor

mid_ps$salinity <- as.factor(mid_ps$salinity)

Assign treatment as characters from integers then to factors

mid_ps$nitrate <- as.factor(as.character(mid_ps$treatment))

Assign deltaNPQ as a numeric

mid_ps$deltaNPQ <- as.numeric(mid_ps$deltaNPQ)

Assign new column for chronological ranking of individuals in each day_group

mid_ps <- mid_ps %>%
  group_by(day_group) %>%
  mutate(rank = rank(rlc_end_time)) %>%
  ungroup()

Use ranking to make smaller groups of ~15 minutes for rlc_order

mid_ps <- mid_ps %>%
  group_by(day_group) %>%
  mutate(rlc_order = floor((rank -1)/3) +1)

Combine nitrate and salinity for a treatment number

mid_ps <- mid_ps %>%
  mutate(treatment = case_when(
    nitrate == 1 & salinity == 35 ~ "1",
    nitrate == 1 & salinity == 28 ~ "2",
    nitrate == 2 & salinity == 35 ~ "3",
    nitrate == 2 & salinity == 28 ~ "4",
    nitrate == 3 & salinity == 35 ~ "5",
    nitrate == 3 & salinity == 28 ~ "6",
    nitrate == 4 & salinity == 35 ~ "7",
    nitrate == 4 & salinity == 28 ~ "8",
  ))

mid_ps_d9 <- subset(mid_ps, rlc_day == 9) #do not separate by plant_part

Toggle between the plant_part for output. Use Day 9 for final analysis

canopy <- subset(mid_ps, rlc_day == 9 & plant_part == "canopy")
canopy$treatment_graph[canopy$treatment == 1] <- "1) 0.5μmol/35 ppt"
canopy$treatment_graph[canopy$treatment == 2] <- "2) 0.5μmol/28 ppt"
canopy$treatment_graph[canopy$treatment == 3] <- "3) 2μmol/35 ppt" 
canopy$treatment_graph[canopy$treatment == 4] <- "4) 2μmol/28 ppt"
canopy$treatment_graph[canopy$treatment == 5] <- "5) 4μmol/35 ppt"
canopy$treatment_graph[canopy$treatment == 6] <- "6) 4μmol/28 ppt"
canopy$treatment_graph[canopy$treatment == 7] <- "7) 8μmol/35 ppt"
canopy$treatment_graph[canopy$treatment == 8] <- "8) 8μmol/28 ppt"

under <- subset(mid_ps, rlc_day == 9 & plant_part == "under")
under$treatment_graph[under$treatment == 1] <- "1) 0.5μmol/35 ppt"
under$treatment_graph[under$treatment == 2] <- "2) 0.5μmol/28 ppt"
under$treatment_graph[under$treatment == 3] <- "3) 2μmol/35 ppt" 
under$treatment_graph[under$treatment == 4] <- "4) 2μmol/28 ppt"
under$treatment_graph[under$treatment == 5] <- "5) 4μmol/35 ppt"
under$treatment_graph[under$treatment == 6] <- "6) 4μmol/28 ppt"
under$treatment_graph[under$treatment == 7] <- "7) 8μmol/35 ppt"
under$treatment_graph[under$treatment == 8] <- "8) 8μmol/28 ppt"

Add new column to subsets from time over 28 summary dataset

time_over_28 <- read_csv("../../../data/midway_2024/transformed/mins_28_plant_part.csv") #load dataset

mean_mins28 <- time_over_28 %>%
  select(plant_part, run, mean_mins28_plant_part) %>% #keep only relevant columns of data
  mutate(run = as.factor(run), mean_mins28_plant_part = as.factor(mean_mins28_plant_part))
  
canopy <- canopy %>%
  left_join(mean_mins28, by = c("run", "plant_part")) #join the datasets

canopy <- canopy %>%
  rename(mean_mins28 = mean_mins28_plant_part) #name is too long

under <- under %>%
  left_join(mean_mins28, by = c("run", "plant_part")) #join the datasets

under <- under %>%
  rename(mean_mins28 = mean_mins28_plant_part) #name is too long

Add growth rate from growth dataset

growth_rate <- read.csv("../../../data/midway_2024/input/midway_growth.csv")
growth_rate$treatment <- as.factor(growth_rate$treatment)

Make a new column for weight change (difference final from initial)

growth_rate$d9_growth_percent <- (growth_rate$final_weight - growth_rate$initial_weight) / growth_rate$initial_weight * 100
d9 <- growth_rate %>%
  select(ID, d9_growth_percent)#keep only relevant columns of data

canopy <- canopy %>%
  left_join(d9, by = c("ID")) #join the datasets   

under <- under %>%
  left_join(d9, by = c("ID"))

Datasets are ready, now let’s run models

Function to run all of the checks with the individual models

check_model_fit <- function(model, terms) {
  print(r.squaredGLMM(model))
  hist(resid(model))
  plot(resid(model) ~ fitted(model))
  qqnorm(resid(model))
  qqline(resid(model))
  print(performance::check_model(model))
  #print(tab_model(model, show.intercept = TRUE, show.se = TRUE, show.stat = TRUE, 
                  #show.df = TRUE, show.zeroinf = TRUE))
  plot(ggpredict(model, terms = terms))
  
}

Get means for each predictor

get_data_means <- function(data, predictors_list, response) {
  data %>%
  group_by(across(all_of(c("nitrate", setdiff(predictors_list, "nitrate"))))) %>%
  summarise(across(all_of(response), list(mean = mean), .group = "drop"))
  
}

Function to run the the models for likelihood ratio tests where REML must be FALSE.

Dataset, response variables, fixed effects (predictors_list), and random effects as arguments.

compare_lmer_models <- function(data, response, predictors_list, random_effects) {
  
  response <- as.character(response) # Ensure response is a string
  
  models <- list() # Create a list to store the models
  
  anova_result <- list() #initialize anova
  
  random_effects_part <- paste0("(1 | ", random_effects, ")", collapse = " + ") #random effects used for all
  
  fixed_effects_all <- paste(predictors_list, collapse = " + ") #make a formula for all predictors
  formula <- as.formula(paste(response, "~", fixed_effects_all, "+", random_effects_part))
  
  model_all <- lmer(formula = formula, data = data, REML = FALSE) #model_all is the version with both predictors for comparison
  
  # Loop through predictors_list to create and fit models
  for (i in seq_along(predictors_list)) {
    fixed_effects <- paste(predictors_list[[i]], collapse = " + ")
    formula <- as.formula(paste(response, "~", fixed_effects, "+", random_effects_part))
    
    models[[i]] <- lmer(formula = formula, data = data, REML = FALSE) # Fit the lmer model and store in the list
    
    anova_result[[i]] <- anova(models[[i]], model_all)
  }
  
  data_means <- get_data_means(data, predictors_list, response)
  
  # Return results as a list
  return(list(
    models = models,
    model_all = model_all,
    anova_result = anova_result,
    data_means = data_means
  ))
}

predictors_list <- c("salinity", "nitrate") #set the predictors list

Run the models with the following inputs

Pmax_____________________________________

Inputs for Canopy/Pmax

canopy_pmax_results <- compare_lmer_models(
  data = canopy,
  response = "pmax",
  predictors_list,
  random_effects = c("plantID", "rlc_order")
)

Access results

summary(canopy_pmax_results$models[[1]]) # Model 1 summary (salinity)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    919.7    932.5   -454.8    909.7       91 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.79668 -0.55604 -0.09338  0.55720  2.60374 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  plantID   (Intercept) 199.9    14.14   
##  rlc_order (Intercept) 157.1    12.53   
##  Residual              462.6    21.51   
## Number of obs: 96, groups:  plantID, 48; rlc_order, 32
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  101.765      4.821  50.311  21.107   <2e-16 ***
## salinity35    -7.549      6.116  39.693  -1.234    0.224    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## salinity35 -0.634

summary(canopy_pmax_results$models[[2]]) # Model 2 summary (nitrate)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    920.2    938.2   -453.1    906.2       89 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.71404 -0.57772 -0.06531  0.66096  2.58147 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  plantID   (Intercept) 169.9    13.03   
##  rlc_order (Intercept) 147.3    12.14   
##  Residual              466.1    21.59   
## Number of obs: 96, groups:  plantID, 48; rlc_order, 32
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  109.257      6.530  50.666  16.732   <2e-16 ***
## nitrate2     -15.133      8.761  47.771  -1.727   0.0906 .  
## nitrate3     -19.788      9.144  50.888  -2.164   0.0352 *  
## nitrate4     -10.144      8.761  47.771  -1.158   0.2527    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) nitrt2 nitrt3
## nitrate2 -0.671              
## nitrate3 -0.700  0.522       
## nitrate4 -0.671  0.455  0.522

summary(canopy_pmax_results$model_all) # Model 3 summary (salinity and nitrate)

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    920.6    941.1   -452.3    904.6       88 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.86204 -0.55814 -0.01586  0.65726  2.69009 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  plantID   (Intercept) 151.4    12.31   
##  rlc_order (Intercept) 154.7    12.44   
##  Residual              463.5    21.53   
## Number of obs: 96, groups:  plantID, 48; rlc_order, 32
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  113.054      7.048  51.337  16.040   <2e-16 ***
## salinity35    -7.581      5.773  39.262  -1.313   0.1967    
## nitrate2     -15.145      8.586  47.480  -1.764   0.0842 .  
## nitrate3     -19.850      8.989  51.181  -2.208   0.0317 *  
## nitrate4     -10.095      8.586  47.480  -1.176   0.2456    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.410                     
## nitrate2   -0.609  0.000              
## nitrate3   -0.638  0.000  0.523       
## nitrate4   -0.609  0.000  0.452  0.523

canopy_pmax_results$anova_result[[1]]        # ANOVA did nitrate have an effect on Pmax?

## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
##             npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## models[[i]]    5 919.70 932.52 -454.85   909.70                     
## model_all      8 920.57 941.09 -452.29   904.57 5.1236  3      0.163

canopy_pmax_results$anova_result[[2]]         # ANOVA did salinity have an effect on Pmax?

## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
##             npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## models[[i]]    7 920.25 938.20 -453.12   906.25                     
## model_all      8 920.57 941.09 -452.29   904.57 1.6766  1     0.1954

check_model_fit(canopy_pmax_results$model_all, terms = predictors_list)

##             R2m     R2c
## [1,] 0.08265697 0.44759

canopy_pmax_results$data_means

## # A tibble: 8 × 3
## # Groups:   nitrate [4]
##   nitrate salinity pmax_mean
##   <fct>   <fct>        <dbl>
## 1 1       28           108. 
## 2 1       35           110. 
## 3 2       28           101. 
## 4 2       35            87.1
## 5 3       28            94.9
## 6 3       35            85.6
## 7 4       28           102. 
## 8 4       35            95.0

Inputs for Understory/Pmax

under_pmax_results <- compare_lmer_models(
  data = under,
  response = "pmax",
  predictors_list,
  random_effects = c("mean_mins28", "rlc_order")
)

Access results

summary(under_pmax_results$models[[1]]) # Model 1 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    883.6    896.5   -436.8    873.6       91 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.39721 -0.58387 -0.02647  0.59074  2.70764 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  rlc_order   (Intercept)  15.89    3.986  
##  mean_mins28 (Intercept)  13.23    3.637  
##  Residual                501.46   22.393  
## Number of obs: 96, groups:  rlc_order, 16; mean_mins28, 2
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept) 75.62709    4.26002  4.33413  17.753 3.26e-05 ***
## salinity35   0.05666    4.61214 90.54020   0.012     0.99    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## salinity35 -0.541

summary(under_pmax_results$models[[2]]) # Model 2 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    886.8    904.7   -436.4    872.8       89 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.37932 -0.63396 -0.05697  0.61647  2.84997 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  rlc_order   (Intercept)  12.06    3.473  
##  mean_mins28 (Intercept)  13.05    3.613  
##  Residual                500.01   22.361  
## Number of obs: 96, groups:  rlc_order, 16; mean_mins28, 2
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)   79.057      5.383  9.864  14.687 4.99e-08 ***
## nitrate2      -6.098      6.598 79.586  -0.924    0.358    
## nitrate3      -4.421      6.694 45.972  -0.660    0.512    
## nitrate4      -3.088      6.598 79.586  -0.468    0.641    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) nitrt2 nitrt3
## nitrate2 -0.613              
## nitrate3 -0.622  0.507       
## nitrate4 -0.613  0.485  0.507

summary(under_pmax_results$model_all) # Model 3 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    888.8    909.3   -436.4    872.8       88 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.37809 -0.63522 -0.05697  0.61709  2.85120 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  rlc_order   (Intercept)  12.06    3.472  
##  mean_mins28 (Intercept)  13.05    3.613  
##  Residual                500.01   22.361  
## Number of obs: 96, groups:  rlc_order, 16; mean_mins28, 2
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept) 79.02837    5.85285 13.55302  13.503 3.03e-09 ***
## salinity35   0.05736    4.59678 90.32593   0.012    0.990    
## nitrate2    -6.09755    6.59829 79.58591  -0.924    0.358    
## nitrate3    -4.42113    6.69434 45.97127  -0.660    0.512    
## nitrate4    -3.08787    6.59829 79.58591  -0.468    0.641    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.393                     
## nitrate2   -0.564  0.000              
## nitrate3   -0.572  0.000  0.507       
## nitrate4   -0.564  0.000  0.485  0.507

under_pmax_results$anova_result[[1]]        # ANOVA did nitrate have an effect on Pmax?

## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
##             npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## models[[i]]    5 883.65 896.47 -436.82   873.65                     
## model_all      8 888.75 909.27 -436.38   872.75 0.8926  3     0.8272

under_pmax_results$anova_result[[2]]         # ANOVA did salinity have an effect on Pmax?

## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
##             npar    AIC    BIC  logLik deviance Chisq Df Pr(>Chisq)
## models[[i]]    7 886.75 904.70 -436.38   872.75                    
## model_all      8 888.75 909.27 -436.38   872.75 2e-04  1       0.99

check_model_fit(under_pmax_results$model_all, terms = predictors_list)

##              R2m        R2c
## [1,] 0.009520906 0.05687984

under_pmax_results$data_means

## # A tibble: 8 × 3
## # Groups:   nitrate [4]
##   nitrate salinity pmax_mean
##   <fct>   <fct>        <dbl>
## 1 1       28            78.1
## 2 1       35            80.4
## 3 2       28            73.0
## 4 2       35            73.5
## 5 3       28            68.7
## 6 3       35            79.6
## 7 4       28            82.7
## 8 4       35            69.2

NPQmax______________________________

Inputs for Canopy/NPQmax

canopy_npqmax_results <- compare_lmer_models(
  data = canopy,
  response = "maxNPQ_Ypoint1",
  predictors_list,
  random_effects = c("plantID", "mean_mins28", "rlc_order")
)

Access results

summary(canopy_npqmax_results$models[[1]]) # Model 1 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##     65.0     80.4    -26.5     53.0       90 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5570 -0.4536 -0.1271  0.4531  4.0118 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  plantID     (Intercept) 0.042708 0.20666 
##  rlc_order   (Intercept) 0.001894 0.04352 
##  mean_mins28 (Intercept) 0.012323 0.11101 
##  Residual                0.062045 0.24909 
## Number of obs: 96, groups:  plantID, 48; rlc_order, 32; mean_mins28, 2
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)   
## (Intercept)  0.71940    0.09643  3.59648   7.460   0.0026 **
## salinity35  -0.03194    0.07855 46.01523  -0.407   0.6862   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## salinity35 -0.407

summary(canopy_npqmax_results$models[[2]]) # Model 2 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##     42.4     62.9    -13.2     26.4       88 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1222 -0.4656  0.0082  0.4043  5.2251 
## 
## Random effects:
##  Groups      Name        Variance  Std.Dev.
##  plantID     (Intercept) 0.0105110 0.10252 
##  rlc_order   (Intercept) 0.0008604 0.02933 
##  mean_mins28 (Intercept) 0.0117902 0.10858 
##  Residual                0.0630749 0.25115 
## Number of obs: 96, groups:  plantID, 48; rlc_order, 32; mean_mins28, 2
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  1.00567    0.09725  4.06386  10.341 0.000455 ***
## nitrate2    -0.36564    0.08413 46.87709  -4.346 7.40e-05 ***
## nitrate3    -0.37748    0.08442 41.65616  -4.472 5.87e-05 ***
## nitrate4    -0.46583    0.08413 46.87709  -5.537 1.35e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) nitrt2 nitrt3
## nitrate2 -0.433              
## nitrate3 -0.434  0.502       
## nitrate4 -0.433  0.497  0.502

summary(canopy_npqmax_results$model_all) # Model 3 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##     44.1     67.1    -13.0     26.1       87 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1856 -0.4654  0.0101  0.3926  5.1965 
## 
## Random effects:
##  Groups      Name        Variance  Std.Dev.
##  plantID     (Intercept) 0.0102212 0.1011  
##  rlc_order   (Intercept) 0.0006758 0.0260  
##  mean_mins28 (Intercept) 0.0117925 0.1086  
##  Residual                0.0632712 0.2515  
## Number of obs: 96, groups:  plantID, 48; rlc_order, 32; mean_mins28, 2
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  1.02164    0.10151  4.79409  10.064  0.00021 ***
## salinity35  -0.03210    0.05914 45.32554  -0.543  0.58987    
## nitrate2    -0.36563    0.08386 46.79540  -4.360 7.08e-05 ***
## nitrate3    -0.37723    0.08408 41.28225  -4.487 5.68e-05 ***
## nitrate4    -0.46578    0.08386 46.79540  -5.555 1.28e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.291                     
## nitrate2   -0.413  0.000              
## nitrate3   -0.414  0.000  0.501       
## nitrate4   -0.413  0.000  0.497  0.501

canopy_npqmax_results$anova[[1]]        # ANOVA did nitrate have an effect on NPQmax?

## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
##             npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## models[[i]]    6 64.984 80.371 -26.492   52.984                         
## model_all      9 44.068 67.147 -13.034   26.068 26.916  3  6.131e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

canopy_npqmax_results$anova[[2]]        # ANOVA did salinity have an effect on NPQmax?

## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
##             npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## models[[i]]    8 42.362 62.876 -13.181   26.362                     
## model_all      9 44.068 67.147 -13.034   26.068 0.2933  1     0.5881

check_model_fit(canopy_npqmax_results$model_all, terms = predictors_list)

##            R2m       R2c
## [1,] 0.2745407 0.4660269

canopy_npqmax_results$data_means

## # A tibble: 8 × 3
## # Groups:   nitrate [4]
##   nitrate salinity maxNPQ_Ypoint1_mean
##   <fct>   <fct>                  <dbl>
## 1 1       28                     1.01 
## 2 1       35                     1.00 
## 3 2       28                     0.619
## 4 2       35                     0.660
## 5 3       28                     0.704
## 6 3       35                     0.554
## 7 4       28                     0.546
## 8 4       35                     0.534

Inputs for Understory/NPQmax

under_npqmax_results <- compare_lmer_models(
  data = under,
  response = "maxNPQ_Ypoint1",
  predictors_list,
  random_effects = c("plantID")
)

Access results

summary(under_npqmax_results$models[[1]]) # Model 1 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    -17.1     -6.8     12.5    -25.1       92 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.68725 -0.49727  0.06191  0.44741  2.92434 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plantID  (Intercept) 0.00591  0.07688 
##  Residual             0.03956  0.19889 
## Number of obs: 96, groups:  plantID, 48
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  0.57165    0.03272 48.00000  17.473   <2e-16 ***
## salinity35   0.08319    0.04627 48.00000   1.798   0.0785 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## salinity35 -0.707

summary(under_npqmax_results$models[[2]]) # Model 2 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    -35.3    -19.9     23.7    -47.3       90 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.09356 -0.48225  0.08922  0.55435  2.49521 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plantID  (Intercept) 0.00000  0.0000  
##  Residual             0.03577  0.1891  
## Number of obs: 96, groups:  plantID, 48
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  0.79508    0.03861 96.00000  20.595  < 2e-16 ***
## nitrate2    -0.25063    0.05460 96.00000  -4.590 1.34e-05 ***
## nitrate3    -0.21083    0.05460 96.00000  -3.862 0.000204 ***
## nitrate4    -0.26592    0.05460 96.00000  -4.871 4.38e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) nitrt2 nitrt3
## nitrate2 -0.707              
## nitrate3 -0.707  0.500       
## nitrate4 -0.707  0.500  0.500
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

summary(under_npqmax_results$model_all) # Model 3 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    -38.1    -20.1     26.0    -52.1       89 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0301 -0.5251  0.1471  0.6169  2.7833 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plantID  (Intercept) 0.00000  0.0000  
##  Residual             0.03404  0.1845  
## Number of obs: 96, groups:  plantID, 48
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  0.75349    0.04211 96.00000  17.895  < 2e-16 ***
## salinity35   0.08319    0.03766 96.00000   2.209 0.029563 *  
## nitrate2    -0.25063    0.05326 96.00000  -4.706 8.50e-06 ***
## nitrate3    -0.21083    0.05326 96.00000  -3.959 0.000145 ***
## nitrate4    -0.26592    0.05326 96.00000  -4.993 2.66e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.447                     
## nitrate2   -0.632  0.000              
## nitrate3   -0.632  0.000  0.500       
## nitrate4   -0.632  0.000  0.500  0.500
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

under_npqmax_results$anova[[1]]         # ANOVA did nitrate have an effect on NPQmax?

## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
##             npar     AIC      BIC logLik deviance  Chisq Df Pr(>Chisq)    
## models[[i]]    4 -17.097  -6.8392 12.548  -25.097                         
## model_all      7 -38.065 -20.1147 26.033  -52.065 26.969  3  5.977e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

under_npqmax_results$anova[[2]]         # ANOVA did salinity have an effect on NPQmax?

## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
##             npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)  
## models[[i]]    6 -35.306 -19.920 23.653  -47.306                       
## model_all      7 -38.065 -20.115 26.033  -52.065 4.7592  1    0.02914 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_model_fit(under_npqmax_results$model_all, terms = predictors_list)

##            R2m       R2c
## [1,] 0.2808741 0.2808741

under_npqmax_results$data_means

## # A tibble: 8 × 3
## # Groups:   nitrate [4]
##   nitrate salinity maxNPQ_Ypoint1_mean
##   <fct>   <fct>                  <dbl>
## 1 1       28                     0.724
## 2 1       35                     0.866
## 3 2       28                     0.509
## 4 2       35                     0.580
## 5 3       28                     0.550
## 6 3       35                     0.619
## 7 4       28                     0.504
## 8 4       35                     0.554

deltaNPQ____________________________

Inputs for Canopy/deltaNPQ

canopy_delta_npq_results <- compare_lmer_models(
  data = canopy,
  response = "deltaNPQ",
  predictors_list,
  random_effects = c("plantID")
)

Access results

summary(canopy_delta_npq_results$models[[1]]) # Model 1 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##     47.9     58.1    -19.9     39.9       92 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6774 -0.5045 -0.1333  0.3854  4.1555 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plantID  (Intercept) 0.03378  0.1838  
##  Residual             0.06112  0.2472  
## Number of obs: 96, groups:  plantID, 48
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  0.53806    0.05178 48.00000  10.392 7.09e-14 ***
## salinity35  -0.02440    0.07322 48.00000  -0.333     0.74    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## salinity35 -0.707

summary(canopy_delta_npq_results$models[[2]]) # Model 2 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##     24.4     39.8     -6.2     12.4       90 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6222 -0.5588 -0.1118  0.4499  5.3800 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plantID  (Intercept) 0.005754 0.07585 
##  Residual             0.061116 0.24722 
## Number of obs: 96, groups:  plantID, 48
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  0.81296    0.05501 48.00000  14.779  < 2e-16 ***
## nitrate2    -0.36958    0.07779 48.00000  -4.751 1.87e-05 ***
## nitrate3    -0.35367    0.07779 48.00000  -4.546 3.71e-05 ***
## nitrate4    -0.42513    0.07779 48.00000  -5.465 1.63e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) nitrt2 nitrt3
## nitrate2 -0.707              
## nitrate3 -0.707  0.500       
## nitrate4 -0.707  0.500  0.500

summary(canopy_delta_npq_results$model_all) # Model 3 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##     26.2     44.2     -6.1     12.2       89 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5826 -0.5452 -0.0900  0.4251  5.3485 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plantID  (Intercept) 0.005605 0.07487 
##  Residual             0.061116 0.24722 
## Number of obs: 96, groups:  plantID, 48
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  0.82516    0.06138 48.00000  13.444  < 2e-16 ***
## salinity35  -0.02440    0.05490 48.00000  -0.444    0.659    
## nitrate2    -0.36958    0.07763 48.00000  -4.761 1.81e-05 ***
## nitrate3    -0.35367    0.07763 48.00000  -4.556 3.60e-05 ***
## nitrate4    -0.42513    0.07763 48.00000  -5.476 1.57e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.447                     
## nitrate2   -0.632  0.000              
## nitrate3   -0.632  0.000  0.500       
## nitrate4   -0.632  0.000  0.500  0.500

canopy_delta_npq_results$anova[[1]]         # ANOVA did nitrate have an effect on deltaNPQ?

## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
##             npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## models[[i]]    4 47.858 58.115 -19.929   39.858                         
## model_all      7 26.202 44.152  -6.101   12.202 27.656  3   4.29e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

canopy_delta_npq_results$anova[[2]]         # ANOVA did salinity have an effect on deltaNPQ?

## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
##             npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## models[[i]]    6 24.399 39.785 -6.1995   12.399                     
## model_all      7 26.202 44.152 -6.1010   12.202 0.1971  1     0.6571

check_model_fit(canopy_delta_npq_results$model_all, terms = predictors_list)

##            R2m       R2c
## [1,] 0.3002212 0.3590055

canopy_delta_npq_results$data_means

## # A tibble: 8 × 3
## # Groups:   nitrate [4]
##   nitrate salinity deltaNPQ_mean
##   <fct>   <fct>            <dbl>
## 1 1       28               0.814
## 2 1       35               0.812
## 3 2       28               0.443
## 4 2       35               0.444
## 5 3       28               0.514
## 6 3       35               0.405
## 7 4       28               0.382
## 8 4       35               0.394

Inputs for Understory/deltaNPQ

under_delta_npq_results <- compare_lmer_models(
  data = under,
  response = "deltaNPQ",
  predictors_list,
  random_effects = c("plantID")
)

Access results

summary(under_delta_npq_results$models[[1]]) # Model 1 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    -59.4    -49.1     33.7    -67.4       92 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.15735 -0.57743 -0.01169  0.42563  2.39882 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plantID  (Intercept) 0.01097  0.1047  
##  Residual             0.02006  0.1416  
## Number of obs: 96, groups:  plantID, 48
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  0.39404    0.02958 48.00000   13.32   <2e-16 ***
## salinity35   0.06402    0.04183 48.00000    1.53    0.132    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## salinity35 -0.707

summary(under_delta_npq_results$models[[2]]) # Model 2 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    -95.1    -79.7     53.5   -107.1       90 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1095 -0.4721  0.1215  0.4697  2.2747 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plantID  (Intercept) 0.0000   0.0000  
##  Residual             0.0192   0.1386  
## Number of obs: 96, groups:  plantID, 48
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  0.62083    0.02828 96.00000  21.952  < 2e-16 ***
## nitrate2    -0.26313    0.04000 96.00000  -6.579 2.49e-09 ***
## nitrate3    -0.23767    0.04000 96.00000  -5.942 4.51e-08 ***
## nitrate4    -0.27833    0.04000 96.00000  -6.959 4.21e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) nitrt2 nitrt3
## nitrate2 -0.707              
## nitrate3 -0.707  0.500       
## nitrate4 -0.707  0.500  0.500
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

summary(under_delta_npq_results$model_all) # Model 3 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    -98.3    -80.4     56.2   -112.3       89 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.95853 -0.39173  0.03624  0.57564  2.10050 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  plantID  (Intercept) 0.00000  0.0000  
##  Residual             0.01817  0.1348  
## Number of obs: 96, groups:  plantID, 48
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  0.58882    0.03076 96.00000  19.140  < 2e-16 ***
## salinity35   0.06402    0.02752 96.00000   2.327   0.0221 *  
## nitrate2    -0.26313    0.03891 96.00000  -6.762 1.06e-09 ***
## nitrate3    -0.23767    0.03891 96.00000  -6.107 2.15e-08 ***
## nitrate4    -0.27833    0.03891 96.00000  -7.152 1.69e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.447                     
## nitrate2   -0.632  0.000              
## nitrate3   -0.632  0.000  0.500       
## nitrate4   -0.632  0.000  0.500  0.500
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

under_delta_npq_results$anova[[1]]         # ANOVA did nitrate have an effect on deltaNPQ?

## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
##             npar     AIC     BIC logLik deviance Chisq Df Pr(>Chisq)    
## models[[i]]    4 -59.359 -49.102 33.679  -67.359                        
## model_all      7 -98.319 -80.368 56.159 -112.319 44.96  3  9.437e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

under_delta_npq_results$anova[[2]]         # ANOVA did salinity have an effect on deltaNPQ?

## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
##             npar     AIC     BIC logLik deviance Chisq Df Pr(>Chisq)  
## models[[i]]    6 -95.053 -79.667 53.526  -107.05                      
## model_all      7 -98.319 -80.368 56.159  -112.32 5.266  1    0.02175 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

check_model_fit(under_delta_npq_results$model_all, terms = predictors_list)

##           R2m      R2c
## [1,] 0.435658 0.435658

under_delta_npq_results$data_means

## # A tibble: 8 × 3
## # Groups:   nitrate [4]
##   nitrate salinity deltaNPQ_mean
##   <fct>   <fct>            <dbl>
## 1 1       28               0.543
## 2 1       35               0.699
## 3 2       28               0.351
## 4 2       35               0.365
## 5 3       28               0.350
## 6 3       35               0.417
## 7 4       28               0.333
## 8 4       35               0.352

Ek____________________________________________

Inputs for Canopy/Ek

canopy_ek_results <- compare_lmer_models(
  data = canopy,
  response = "ek.est",
  predictors_list,
  random_effects = c("plantID", "mean_mins28", "rlc_order")
)

Access results

summary(canopy_ek_results$models[[1]]) # Model 1 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    993.6   1009.0   -490.8    981.6       90 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.8311 -0.5932 -0.0605  0.5266  3.3316 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  plantID     (Intercept)  361.2   19.01   
##  rlc_order   (Intercept)  240.9   15.52   
##  mean_mins28 (Intercept)  150.7   12.28   
##  Residual                1040.1   32.25   
## Number of obs: 96, groups:  plantID, 48; rlc_order, 32; mean_mins28, 2
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  177.380     10.964   3.239  16.179 0.000329 ***
## salinity35   -13.284      8.712  39.775  -1.525 0.135214    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## salinity35 -0.397

summary(canopy_ek_results$models[[2]]) # Model 2 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    995.7   1016.2   -489.8    979.7       88 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.1327 -0.5942 -0.0610  0.4520  3.0012 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  plantID     (Intercept)  316.0   17.78   
##  rlc_order   (Intercept)  258.3   16.07   
##  mean_mins28 (Intercept)  148.9   12.20   
##  Residual                1033.6   32.15   
## Number of obs: 96, groups:  plantID, 48; rlc_order, 32; mean_mins28, 2
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  159.526     12.649   5.654  12.612 2.36e-05 ***
## nitrate2      12.348     12.498  46.880   0.988   0.3282    
## nitrate3       7.310     12.981  50.099   0.563   0.5759    
## nitrate4      25.190     12.498  46.880   2.016   0.0496 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) nitrt2 nitrt3
## nitrate2 -0.494              
## nitrate3 -0.513  0.519       
## nitrate4 -0.494  0.461  0.519

summary(canopy_ek_results$model_all) # Model 3 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    995.0   1018.1   -488.5    977.0       87 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.94436 -0.66754  0.00851  0.56066  3.12130 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  plantID     (Intercept)  247.8   15.74   
##  rlc_order   (Intercept)  301.7   17.37   
##  mean_mins28 (Intercept)  145.2   12.05   
##  Residual                1020.1   31.94   
## Number of obs: 96, groups:  plantID, 48; rlc_order, 32; mean_mins28, 2
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  166.014     13.063   6.446  12.709 8.39e-06 ***
## salinity35   -13.615      8.118  36.394  -1.677   0.1021    
## nitrate2      12.618     12.067  45.994   1.046   0.3012    
## nitrate3       7.773     12.626  50.791   0.616   0.5409    
## nitrate4      25.734     12.067  45.994   2.133   0.0383 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.311                     
## nitrate2   -0.462  0.000              
## nitrate3   -0.483  0.000  0.523       
## nitrate4   -0.462  0.000  0.453  0.523

canopy_ek_results$anova[[1]]         # ANOVA did nitrate have an effect on Ek?

## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
##             npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## models[[i]]    6 993.60 1009.0 -490.80   981.60                     
## model_all      9 995.04 1018.1 -488.52   977.04 4.5637  3     0.2067

canopy_ek_results$anova[[2]]         # ANOVA did salinity have an effect on Ek?

## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
##             npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## models[[i]]    8 995.68 1016.2 -489.84   979.68                     
## model_all      9 995.04 1018.1 -488.52   977.04 2.6442  1     0.1039

check_model_fit(canopy_ek_results$model_all, terms = predictors_list)

##             R2m      R2c
## [1,] 0.07310292 0.448615

canopy_ek_results$data_means

## # A tibble: 8 × 3
## # Groups:   nitrate [4]
##   nitrate salinity ek.est_mean
##   <fct>   <fct>          <dbl>
## 1 1       28              158.
## 2 1       35              164.
## 3 2       28              178.
## 4 2       35              166.
## 5 3       28              176.
## 6 3       35              157.
## 7 4       28              196.
## 8 4       35              171.

Inputs for Understory/Ek

under_ek_results <- compare_lmer_models(
  data = under,
  response = "ek.est",
  predictors_list,
  random_effects = c("mean_mins28", "rlc_order")
)

Access results

summary(under_ek_results$models[[1]]) # Model 1 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    970.8    983.7   -480.4    960.8       91 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9867 -0.6906 -0.1270  0.6353  3.0333 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  rlc_order   (Intercept)   36.832  6.069  
##  mean_mins28 (Intercept)    2.534  1.592  
##  Residual                1264.120 35.554  
## Number of obs: 96, groups:  rlc_order, 16; mean_mins28, 2
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  133.936      5.489   6.557  24.402  1.1e-07 ***
## salinity35    -6.767      7.318  89.907  -0.925    0.358    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## salinity35 -0.667

summary(under_ek_results$models[[2]]) # Model 2 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    969.7    987.7   -477.9    955.7       89 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.12670 -0.68008 -0.08205  0.63288  2.99702 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  rlc_order   (Intercept)   99.844  9.992  
##  mean_mins28 (Intercept)    7.511  2.741  
##  Residual                1144.211 33.826  
## Number of obs: 96, groups:  rlc_order, 16; mean_mins28, 2
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  114.119      7.986  17.418  14.291 4.68e-11 ***
## nitrate2      19.059     10.416  85.624   1.830   0.0708 .  
## nitrate3      25.003     10.856  55.420   2.303   0.0250 *  
## nitrate4      21.671     10.416  85.624   2.081   0.0405 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr) nitrt2 nitrt3
## nitrate2 -0.652              
## nitrate3 -0.680  0.521       
## nitrate4 -0.652  0.457  0.521

summary(under_ek_results$model_all) # Model 3 summary

## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: formula
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    970.8    991.3   -477.4    954.8       88 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.22521 -0.76232 -0.08919  0.57191  2.91970 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  rlc_order   (Intercept)  100.788 10.039  
##  mean_mins28 (Intercept)    7.803  2.793  
##  Residual                1131.695 33.641  
## Number of obs: 96, groups:  rlc_order, 16; mean_mins28, 2
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)  117.457      8.703  23.877  13.496 1.15e-12 ***
## salinity35    -6.748      7.004  86.873  -0.964   0.3379    
## nitrate2      19.080     10.369  85.738   1.840   0.0692 .  
## nitrate3      25.086     10.813  55.664   2.320   0.0240 *  
## nitrate4      21.713     10.369  85.738   2.094   0.0392 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr) slnt35 nitrt2 nitrt3
## salinity35 -0.402                     
## nitrate2   -0.596  0.000              
## nitrate3   -0.621  0.000  0.521       
## nitrate4   -0.596  0.000  0.456  0.521

under_ek_results$anova[[1]]         # ANOVA did nitrate have an effect on Ek?

## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
##             npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## models[[i]]    5 970.83 983.65 -480.41   960.83                     
## model_all      8 970.79 991.30 -477.39   954.79 6.0408  3     0.1096

under_ek_results$anova[[2]]         # ANOVA did salinity have an effect on Ek?

## Data: data
## Models:
## models[[i]]: formula
## model_all: formula
##             npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)
## models[[i]]    7 969.71 987.66 -477.86   955.71                     
## model_all      8 970.79 991.30 -477.39   954.79 0.9237  1     0.3365

check_model_fit(under_ek_results$model_all, terms = predictors_list)

##             R2m       R2c
## [1,] 0.07972884 0.1603012

under_ek_results$data_means

## # A tibble: 8 × 3
## # Groups:   nitrate [4]
##   nitrate salinity ek.est_mean
##   <fct>   <fct>          <dbl>
## 1 1       28              122.
## 2 1       35              111.
## 3 2       28              133.
## 4 2       35              136.
## 5 3       28              133.
## 6 3       35              139.
## 7 4       28              148.
## 8 4       35              123.

ALL OF THE ABOVE RANDOM EFFECTS HAVE BEEN CHECKED AND MAXIMIZED FOR FIT

Plots for raw data

HISTOGRAMS and PLOTS_____________________________________

Function for raw data plots

raw_plots <- function(data, response, label, pretty_color, aescolor, x, y, title, subtitle, ylim1, ylim2, 
                      color1, color2, yint, vjust_t, hjust_t, vjust_s, hjust_s, labels1) {
  histo <- data %>%
    ggplot(aes(x = {{response}})) +
    geom_histogram(fill = pretty_color, color = "black", size = 0.25, alpha = 0.85) +
    labs(title = label) +
    theme_bw()

  plot <- data %>% 
    ggplot(aes(x = treatment_graph, y = {{response}})) + 
    geom_boxplot(size=0.5) + 
    geom_point(alpha = 0.75, size = 3, aes(color = {{aescolor}}), position = "jitter", show.legend = TRUE) + 
    labs(x = x, y = y, title = title, subtitle = subtitle) + 
    #scale_x_discrete(labels = c("0.5μmolN", "2μmolN", "4μmolN", "8μmolN")) + 
    ylim(ylim1, ylim2) + stat_mean() + 
    scale_color_manual(values = c(color1, color2)) +
    scale_x_discrete(labels = labels1) +
    geom_hline(yintercept=yint, color = "red", size = 0.5, alpha = 0.5) +
    theme_bw() +
    theme(legend.position = c(0.90,0.90), plot.title = element_text(face = "bold", vjust = vjust_t, hjust = hjust_t), 
          plot.subtitle = element_text(face = "italic", size = 14, vjust = vjust_s, hjust = hjust_s))
  
  return(list(
    histo = histo,
    plot = plot)) 
}

INPUTS for histograms and ggplots—————–

Pmax ————————

Inputs for canopy/pmax

pmax_canopy <- raw_plots(
  data = canopy,
  response = pmax,
  label = "Canopy",
  pretty_color = "goldenrod1",
  aescolor = salinity,
  x = "Nitrate",
  y = "Day 9 Pmax (μmols electrons m-2 s-1)",
  title = "A - Canopy",
  subtitle = "Chondria tumulosa",
  ylim1 = 0, 
  ylim2 = 175,
  color1 = "goldenrod3",
  color2 = "darkgoldenrod1",
  yint = 100,
  vjust_t = -15,
  hjust_t = 0.05,
  vjust_s = -20,
  hjust_s = 0.05,
  labels1 = c("0.5 μmol", "0.5 μmol", 
              "2 μmol", "2 μmol",
              "4 μmol", "4 μmol",
              "8 μmol", "8 μmol")
)
plot(pmax_canopy$histo)

plot(pmax_canopy$plot)

Inputs for under/pmax

pmax_under <- raw_plots(
  data = under,
  response = pmax,
  label = "Understory",
  pretty_color = "maroon4",
  aescolor = salinity,
  x = "Nitrate",
  y = "Day 9 Pmax (μmols electrons m-2 s-1)",
  title = "B - Understory",
  subtitle = "Chondria tumulosa",
  ylim1 = 0, 
  ylim2 = 175,
  color1 = "deeppink4",
  color2 = "deeppink3",
  yint = 100,
  vjust_t = -15,
  hjust_t = 0.05,
  vjust_s = -20,
  hjust_s = 0.05,
  labels1 = c("0.5 μmol", "0.5 μmol", 
              "2 μmol", "2 μmol",
              "4 μmol", "4 μmol",
              "8 μmol", "8 μmol")
)
plot(pmax_under$histo)

plot(pmax_under$plot)

NPQmax———————

Inputs for canopy/NPQmax

npqmax_canopy <- raw_plots(
  data = canopy,
  response = maxNPQ_Ypoint1,
  label = "Canopy",
  pretty_color = "goldenrod2",
  aescolor = salinity,
  x = "Nitrate",
  y = "Day 9 NPQmax",
  title = "C - Canopy",
  subtitle = "Chondria tumulosa",
  ylim1 = 0, 
  ylim2 = 2,
  color1 = "gold1",
  color2 = "gold3",
  yint = 1,
  vjust_t = -15,
  hjust_t = 0.05,
  vjust_s = -20,
  hjust_s = 0.05,
  labels1 = c("0.5 μmol", "0.5 μmol", 
              "2 μmol", "2 μmol",
              "4 μmol", "4 μmol",
              "8 μmol", "8 μmol")
)
plot(npqmax_canopy$histo)

plot(npqmax_canopy$plot)

Inputs for under/NPQmax

npqmax_under <- raw_plots(
  data = under,
  response = maxNPQ_Ypoint1,
  label = "Understory",
  pretty_color = "maroon1",
  aescolor = salinity,
  x = "Nitrate",
  y = "Day 9 NPQmax",
  title = "D - Understory",
  subtitle = "Chondria tumulosa",
  ylim1 = 0, 
  ylim2 = 2,
  color1 = "hotpink2",
  color2 = "hotpink4",
  yint = 1,
  vjust_t = -15,
  hjust_t = 0.05,
  vjust_s = -20,
  hjust_s = 0.05,
  labels1 = c("0.5 μmol", "0.5 μmol", 
              "2 μmol", "2 μmol",
              "4 μmol", "4 μmol",
              "8 μmol", "8 μmol")
)
plot(npqmax_under$histo)

plot(npqmax_under$plot)

deltaNPQ——————-

Inputs for canopy/deltaNPQ

delta_npq_canopy <- raw_plots(
  data = canopy,
  response = deltaNPQ,
  label = "Canopy",
  pretty_color = "goldenrod3",
  aescolor = salinity,
  x = "Nitrate",
  y = "Day 9 deltaNPQ",
  title = "E - Canopy",
  subtitle = "Chondria tumulosa",
  ylim1 = 0, 
  ylim2 = 2,
  color1 = "lightgoldenrod2",
  color2 = "lightgoldenrod3",
  yint = 0.5,
  vjust_t = -15,
  hjust_t = 0.05,
  vjust_s = -20,
  hjust_s = 0.05,
  labels1 = c("0.5 μmol", "0.5 μmol", 
              "2 μmol", "2 μmol",
              "4 μmol", "4 μmol",
              "8 μmol", "8 μmol")
)
plot(delta_npq_canopy$histo)

plot(delta_npq_canopy$plot)

Inputs for under/deltaNPQ

delta_npq_under <- raw_plots(
  data = under,
  response = deltaNPQ,
  label = "Understory",
  pretty_color = "maroon3",
  aescolor = salinity,
  x = "Nitrate",
  y = "Day 9 deltaNPQ",
  title = "F - Understory",
  subtitle = "Chondria tumulosa",
  ylim1 = 0, 
  ylim2 = 2,
  color1 = "hotpink1",
  color2 = "hotpink3",
  yint = 0.5,
  vjust_t = -15,
  hjust_t = 0.05,
  vjust_s = -20,
  hjust_s = 0.05,
  labels1 = c("0.5 μmol", "0.5 μmol", 
              "2 μmol", "2 μmol",
              "4 μmol", "4 μmol",
              "8 μmol", "8 μmol")
)
plot(delta_npq_under$histo)

plot(delta_npq_under$plot)

Ek————————

Inputs for canopy/Ek

ek_canopy <- raw_plots(
  data = canopy,
  response = ek.est,
  label = "Canopy",
  pretty_color = "goldenrod4",
  aescolor = salinity,
  x = "Nitrate",
  y = "Day 9 Ek (μmols electrons m-2 s-1)",
  title = "G - Canopy",
  subtitle = "Chondria tumulosa",
  ylim1 = 50, 
  ylim2 = 300,
  color1 = "khaki4",
  color2 = "khaki3",
  yint = 150,
  vjust_t = -15,
  hjust_t = 0.05,
  vjust_s = -20,
  hjust_s = 0.05,
  labels1 = c("0.5 μmol", "0.5 μmol", 
              "2 μmol", "2 μmol",
              "4 μmol", "4 μmol",
              "8 μmol", "8 μmol")
)
plot(ek_canopy$histo)

plot(ek_canopy$plot)

Inputs for under/Ek

ek_under <- raw_plots(
  data = under,
  response = ek.est,
  label = "Understory",
  pretty_color = "maroon3",
  aescolor = salinity,
  x = "Nitrate",
  y = "Day 9 Ek (μmols electrons m-2 s-1)",
  title = "H - Understory",
  subtitle = "Chondria tumulosa",
  ylim1 = 50, 
  ylim2 = 300,
  color1 = "palevioletred3",
  color2 = "palevioletred4",
  yint = 150,
  vjust_t = -15,
  hjust_t = 0.05,
  vjust_s = -20,
  hjust_s = 0.05,
  labels1 = c("0.5 μmol", "0.5 μmol", 
              "2 μmol", "2 μmol",
              "4 μmol", "4 μmol",
              "8 μmol", "8 μmol")
)
plot(ek_under$histo)

plot(ek_under$plot)

Midway LMM for Pmax, NPQmax, deltaNPQ, Ek

Angela Richards Donà

2024-12-04

Prepare the data

Load the packages

Load this file for normalized to quantum efficiency of photosynthesis per Silsbe and Kromkamp (2012)

Rename specimenID to simply ID

Standardize date to ISO

Make sure time is hms

Assign run as a factor

Assign temperature as a factor

Assign treatment as characters from integers then to factors

Assign deltaNPQ as a numeric

Assign new column for chronological ranking of individuals in each day_group

Use ranking to make smaller groups of ~15 minutes for rlc_order

Combine nitrate and salinity for a treatment number

Toggle between the plant_part for output. Use Day 9 for final analysis

Add new column to subsets from time over 28 summary dataset

Add growth rate from growth dataset

Make a new column for weight change (difference final from initial)

Datasets are ready, now let’s run models

Function to run all of the checks with the individual models

Get means for each predictor

Function to run the the models for likelihood ratio tests where REML must be FALSE.

Dataset, response variables, fixed effects (predictors_list), and random effects as arguments.

Run the models with the following inputs

Pmax_____________________________________

Inputs for Canopy/Pmax

Access results

Inputs for Understory/Pmax

Access results

NPQmax______________________________

Inputs for Canopy/NPQmax

Access results

Inputs for Understory/NPQmax

Access results

deltaNPQ____________________________

Inputs for Canopy/deltaNPQ

Access results

Inputs for Understory/deltaNPQ

Access results

Ek____________________________________________

Inputs for Canopy/Ek

Access results

Inputs for Understory/Ek

Access results

ALL OF THE ABOVE RANDOM EFFECTS HAVE BEEN CHECKED AND MAXIMIZED FOR FIT

Plots for raw data

HISTOGRAMS and PLOTS_____________________________________

Function for raw data plots

INPUTS for histograms and ggplots—————–

Pmax ————————

Inputs for canopy/pmax

Inputs for under/pmax

NPQmax———————

Inputs for canopy/NPQmax

Inputs for under/NPQmax

deltaNPQ——————-

Inputs for canopy/deltaNPQ

Inputs for under/deltaNPQ

Ek————————

Inputs for canopy/Ek

Inputs for under/Ek